Tarea Evaluación ACP y Cluster

Carga de librerías y datos

Antes de comenzar, se comprueba que las librerías a emplear a lo largo de los ejercicios estén instaladas y, en caso contrario, proceder a instalarlas con el siguiente código:

comprobar <- function(pkg){
    new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
    if (length(new.pkg)) 
        install.packages(new.pkg, dependencies = TRUE)
    sapply(pkg, require, character.only = TRUE)
}
paquetes<-c("tidyverse","factoextra","FactoMineR","cluster","readxl",
            "ggplot2","corrplot","psych","heatmaply","VIM",
            "NbClust","scales","descr","caret","magrittr", "dendextend",
            "ggpubr")
comprobar(paquetes)

##  tidyverse factoextra FactoMineR    cluster     readxl    ggplot2   corrplot 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE 
##      psych  heatmaply        VIM    NbClust     scales      descr      caret 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE 
##   magrittr dendextend     ggpubr 
##       TRUE       TRUE       TRUE

A continuación, antes de proceder a la resolución de las preguntas, se cargarán los datos de las provincias en un dataframe:

# En un script interactivo se podría hacer con el comando choose.files:
# data<-read_xlsx(choose.files(), sheet = "Hoja1")
data<-read_xlsx("Provincias.xlsx", sheet = "Hoja1")
data<-as.data.frame(data)
head(data)

##       Prov Poblacion Mortalidad Natalidad     IPC NumEmpresas Industria
## 1 Albacete    396987       9.41      9.02 101.819       26701      3332
## 2 Alicante   1868438       8.00      8.52 102.309      130438      9992
## 3  Almería    701688       6.91     11.41 101.454       40327      2191
## 4    Álava    321932       7.72     10.26 102.743       19548      2048
## 5 Asturias   1061756      12.16      6.26 102.070       67451      3496
## 6  Badajoz    690929       9.54      8.96 101.036       39640      3059
##   Construccion   CTH Infor  AFS   APT TasaActividad TasaParo Ocupados      PIB
## 1         3428 11111   189  624  3278         57.72    23.28    144.8  7235991
## 2        17418 52720  1851 2745 19865         58.21    23.39    687.1 32139347
## 3         5388 18184   366  944  5840         60.63    31.08    234.4 11909684
## 4         2669  7557   323  354  3341         57.86    13.33    133.0 10716021
## 5         8435 28045   787 1445 10598         51.42    16.97    389.0 21770433
## 6         4599 18420   332  948  5235         56.00    29.63    224.6 10581366
##    CANE     TVF     VS
## 1 20888  215948  30244
## 2 25968 1274096 326705
## 3 22945  395086  72486
## 4  3691  155767   9791
## 5 23910  613905  73250
## 6 37438  372493  49441

Se ha comprobado que los datos cargados son un dataframe por lo que se evaluará su estructura interna para hacerse una idea de su disposición y de los datos faltantes:

summary(data)

##      Prov             Poblacion         Mortalidad       Natalidad     
##  Length:52          Min.   :  84509   Min.   : 5.820   Min.   : 5.550  
##  Class :character   1st Qu.: 322203   1st Qu.: 7.855   1st Qu.: 7.700  
##  Mode  :character   Median : 614723   Median : 9.120   Median : 8.975  
##                     Mean   : 899449   Mean   : 9.379   Mean   : 8.839  
##                     3rd Qu.:1019030   3rd Qu.:10.688   3rd Qu.: 9.610  
##                     Max.   :6454440   Max.   :14.360   Max.   :19.330  
##       IPC         NumEmpresas       Industria      Construccion  
##  Min.   :100.6   Min.   :  3749   Min.   :   75   Min.   :  309  
##  1st Qu.:101.9   1st Qu.: 22822   1st Qu.: 1704   1st Qu.: 2972  
##  Median :102.3   Median : 38000   Median : 2516   Median : 5070  
##  Mean   :102.4   Mean   : 61286   Mean   : 3808   Mean   : 7805  
##  3rd Qu.:102.7   3rd Qu.: 65083   3rd Qu.: 4106   3rd Qu.: 8350  
##  Max.   :104.8   Max.   :508612   Max.   :27416   Max.   :59661  
##       CTH             Infor              AFS               APT        
##  Min.   :  2030   Min.   :   35.0   Min.   :   50.0   Min.   :   504  
##  1st Qu.:  9243   1st Qu.:  185.5   1st Qu.:  486.8   1st Qu.:  3091  
##  Median : 15488   Median :  369.5   Median :  813.5   Median :  5440  
##  Mean   : 23741   Mean   : 1131.9   Mean   : 1378.3   Mean   : 10854  
##  3rd Qu.: 27567   3rd Qu.:  868.2   3rd Qu.: 1429.2   3rd Qu.: 10627  
##  Max.   :158331   Max.   :19058.0   Max.   :12357.0   Max.   :123863  
##  TasaActividad      TasaParo        Ocupados           PIB           
##  Min.   :47.41   Min.   :11.95   Min.   :  24.6   Min.   :  1397441  
##  1st Qu.:55.45   1st Qu.:15.39   1st Qu.: 132.4   1st Qu.:  6509393  
##  Median :57.79   Median :19.51   Median : 222.4   Median : 11883640  
##  Mean   :57.84   Mean   :21.17   Mean   : 347.1   Mean   : 20275145  
##  3rd Qu.:60.07   3rd Qu.:27.75   3rd Qu.: 389.8   3rd Qu.: 21242697  
##  Max.   :68.69   Max.   :37.18   Max.   :2806.4   Max.   :198652445  
##       CANE            TVF                VS        
##  Min.   :    3   Min.   :  26233   Min.   :   200  
##  1st Qu.: 9758   1st Qu.: 211628   1st Qu.: 38941  
##  Median :14037   Median : 335934   Median : 56412  
##  Mean   :19035   Mean   : 484781   Mean   : 70799  
##  3rd Qu.:26020   3rd Qu.: 532066   3rd Qu.: 80625  
##  Max.   :68037   Max.   :2894679   Max.   :326705

aggr_plot <- aggr(data, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(data), cex.axis=.5, gap=1, ylab=c("Histogram of missing data","Pattern"))

## 
##  Variables sorted by number of missings: 
##       Variable Count
##           Prov     0
##      Poblacion     0
##     Mortalidad     0
##      Natalidad     0
##            IPC     0
##    NumEmpresas     0
##      Industria     0
##   Construccion     0
##            CTH     0
##          Infor     0
##            AFS     0
##            APT     0
##  TasaActividad     0
##       TasaParo     0
##       Ocupados     0
##            PIB     0
##           CANE     0
##            TVF     0
##             VS     0

El gráfico muestra que en el dataset no hay datos faltantes, por lo que no será necesario proceder a su inserción.

Por comodidad, se añadirán los nombres de las provincias de la col.1 como rownames eliminando esta columna del conjunto tras la operación.

rownames(data)<-data$Prov
data<-data[,-1]

Pregunta 1: Calcular la matriz de correlaciones, y su representación gráfica ¿Cuáles son las variables más correlacionadas de forma inversa?

Para responder a esta pregunta, se creará la matriz de correlaciones y se representará mediante el paquete corrplot:

matriz_cor<-cor(data, method = "pearson")
corrplot(matriz_cor, type="upper", order="hclust",
         tl.col="black", tl.cex=0.6, tl.srt=90)

Con este gráfico se puede apreciar que las variables más correlacionadas negativamente son Mortalidad con TasaActividad, Natalidad y TasaParo así como la TasaParo con el IPC. Una interpretación de esto es que en las provincias con mayor mortalidad, la natalidad suele ser baja como pudiera ser en las regiones más envejecidas donde hay fallecimientos, pero pocas familias jóvenes.

Pregunta 2: Realizar un análisis de componentes principales sobre la matriz de correlaciones, calculando 7 componentes. Estudiar los valores de los autovalores obtenidos y las gráficas que los resumen. ¿Cuál es el número adecuado de componentes?

Para esta pregunta se utilizará la función PCA de factominer combinada después con fviz de factoextra:

pca<-PCA(data,ncp = 7, scale.unit = TRUE, graph = TRUE)

# Especificando el argumento scale.unit=TRUE el PCA se realiza sobre las variables estandarizadas que es como aplicarlo sobre la matriz de correlaciones

eig<-get_eigenvalue(pca)

Observando los 2 primeros gráficos se extrae que la componente 1 explica el 63,7% de la variabilidad de los datos y la componente 2 el 14,23%, algo muy positivo porque gran parte de la muestra podría quedar reducida a estas 2.

Por otro lado, en cuanto a las provincias como tal, Madrid y Barcelona están muy alejadas del resto de provincias y muy ligadas a la PC1, mientras que hay luego otros 2 grupos más explicados por la PC2, uno con provincias del sur como Melilla, Málaga, Sevilla… y otro con más del centro-norte como Palencia, Pontevedra, Álava, Zamora, etc.

Comprobando los valores de los autovalores:

knitr:: kable(eig, digits = 2, caption="Autovalores")

Autovalores

	eigenvalue	variance.percent	cumulative.variance.percent
Dim.1	11.47	63.70	63.70
Dim.2	2.56	14.23	77.93
Dim.3	1.63	9.08	87.01
Dim.4	0.93	5.19	92.19
Dim.5	0.46	2.54	94.73
Dim.6	0.41	2.30	97.03
Dim.7	0.31	1.71	98.74
Dim.8	0.12	0.65	99.39
Dim.9	0.07	0.41	99.79
Dim.10	0.02	0.11	99.91
Dim.11	0.01	0.05	99.96
Dim.12	0.00	0.02	99.98
Dim.13	0.00	0.01	99.99
Dim.14	0.00	0.00	99.99
Dim.15	0.00	0.00	100.00
Dim.16	0.00	0.00	100.00
Dim.17	0.00	0.00	100.00
Dim.18	0.00	0.00	100.00

Atendiendo a la tabla, con las 4 primeras componentes, se podría explicar el 92,19% de los datos siendo que el valor de la 4 está muy cerca de 1. El número más adecuado de componentes se puede extraer con el scree plot:

fviz_eig(pca, addlabels=TRUE)

El número más óptimo serían 3 o 4 componentes. Como se ha explicado antes, se tomarán 4 para superar el 90% de datos explicados.

Pregunta 3: Hacer de nuevo el análisis sobre la matriz de correlaciones pero ahora indicando el número de componentes principales que hemos decidido retener.

3a: Mostrar los coeficientes para obtener las componentes principales ¿Cuál es la expresión para calcular la primera Componente en función de las variables originales?

En el apartado anterior se concluyó que 4 componentes serían suficientes de manera que:

pca_4<-PCA(data,ncp = 4, scale.unit = TRUE, graph = TRUE)

knitr::kable(pca_4$svd$V, digits=3, caption = "Autovectores")

Autovectores

0.294	0.002	0.050	-0.053
-0.106	-0.527	0.189	-0.161
0.041	0.495	-0.271	-0.110
0.110	-0.365	-0.262	0.435
0.294	-0.026	0.008	-0.069
0.286	-0.045	0.046	0.023
0.293	-0.045	-0.012	-0.026
0.293	-0.011	0.049	-0.028
0.282	-0.042	-0.065	-0.222
0.292	-0.016	0.040	-0.092
0.291	-0.029	-0.028	-0.142
0.114	0.331	-0.363	0.463
-0.014	0.462	0.387	-0.220
0.294	-0.017	0.002	-0.060
0.291	-0.036	-0.037	-0.134
0.018	0.096	0.657	0.278
0.292	-0.002	0.100	0.044
0.172	0.048	0.290	0.567

Cada una de las 4 columnas mostradas representa una de las 4 componenetes escogidas para el análisis y las líneas son las variables originales del dataset con lo que la expresión de la PC1 sería:

$ PC1 = 0,294Poblacion - 0,106Mortalidad + 0,041Natalidad + … + 0,292TVF + 0,172VS $

3b: Mostar una tabla con las correlaciones de las Variables con las Componentes Principales. Para cada Componente indicar las variables con las que está más correlacionada

var<-get_pca_var(pca_4)
knitr::kable(var$cor, digits=2, caption = "Correlaciones PC con variables")

Correlaciones PC con variables

	Dim.1	Dim.2	Dim.3	Dim.4
Poblacion	0.99	0.00	0.06	-0.05
Mortalidad	-0.36	-0.84	0.24	-0.16
Natalidad	0.14	0.79	-0.35	-0.11
IPC	0.37	-0.58	-0.34	0.42
NumEmpresas	1.00	-0.04	0.01	-0.07
Industria	0.97	-0.07	0.06	0.02
Construccion	0.99	-0.07	-0.02	-0.03
CTH	0.99	-0.02	0.06	-0.03
Infor	0.95	-0.07	-0.08	-0.21
AFS	0.99	-0.03	0.05	-0.09
APT	0.98	-0.05	-0.04	-0.14
TasaActividad	0.39	0.53	-0.46	0.45
TasaParo	-0.05	0.74	0.50	-0.21
Ocupados	1.00	-0.03	0.00	-0.06
PIB	0.99	-0.06	-0.05	-0.13
CANE	0.06	0.15	0.84	0.27
TVF	0.99	0.00	0.13	0.04
VS	0.58	0.08	0.37	0.55

De cara a ver las correlaciones más fuertes de las variables con las componentes, se ordenará el df de var$cor de mayor a menor para cada componente y se extraerán las variables con mayor correlación (positiva o negativa):

df<-as.data.frame(var$cor)
df %>% select("Dim.1") %>% arrange(desc(Dim.1))

##                     Dim.1
## Ocupados       0.99699549
## NumEmpresas    0.99616251
## Poblacion      0.99394316
## Construccion   0.99305566
## CTH            0.99170521
## AFS            0.99035020
## TVF            0.98730511
## PIB            0.98507163
## APT            0.98445474
## Industria      0.96717807
## Infor          0.95322382
## VS             0.58359002
## TasaActividad  0.38716977
## IPC            0.37218373
## Natalidad      0.13757149
## CANE           0.06023588
## TasaParo      -0.04739450
## Mortalidad    -0.35992669

En la PC1 tienen una correlación casi perfecta los Ocupados, NumEmpresas, Poblaciónm Construcción, etc. hasta Infor teniendo las correlacionadas negativamente valores muy bajos

df %>% select(Dim.2) %>% arrange(desc(Dim.2))

##                      Dim.2
## Natalidad      0.792749686
## TasaParo       0.738572741
## TasaActividad  0.529035826
## CANE           0.153910432
## VS             0.076521234
## Poblacion      0.003989744
## TVF           -0.002559620
## CTH           -0.016989690
## AFS           -0.026329340
## Ocupados      -0.026944763
## NumEmpresas   -0.041840166
## APT           -0.047097804
## PIB           -0.057858046
## Infor         -0.067235670
## Industria     -0.071587054
## Construccion  -0.072618036
## IPC           -0.584573656
## Mortalidad    -0.843475898

A la PC2 contribuyen positivamente Natalidad, TasaParo y TasaActividad mientras que de forma negativa, Mortalidad e IPC son las que más destacan.

df %>% select(Dim.3) %>% arrange(desc(Dim.3))

##                      Dim.3
## CANE           0.839526520
## TasaParo       0.495165389
## VS             0.370894920
## Mortalidad     0.241432113
## TVF            0.127527357
## Poblacion      0.063848580
## CTH            0.062434517
## Industria      0.059274514
## AFS            0.050562974
## NumEmpresas    0.010165492
## Ocupados       0.003105467
## Construccion  -0.015444759
## APT           -0.035880121
## PIB           -0.047926363
## Infor         -0.082630799
## IPC           -0.335325313
## Natalidad     -0.346016593
## TasaActividad -0.463689989

A la PC3 CANE es la mejor contribuyente en la parte positiva y TasaActividad negativamente aunque sin llegar al 0,5.

df %>% select(Dim.4) %>% arrange(desc(Dim.4))

##                     Dim.4
## VS             0.54810069
## TasaActividad  0.44752178
## IPC            0.42045677
## CANE           0.26905683
## TVF            0.04266413
## Industria      0.02266458
## Construccion  -0.02550892
## CTH           -0.02665589
## Poblacion     -0.05160463
## Ocupados      -0.05803848
## NumEmpresas   -0.06681923
## AFS           -0.08860621
## Natalidad     -0.10599223
## PIB           -0.12912629
## APT           -0.13710328
## Mortalidad    -0.15590654
## TasaParo      -0.21281673
## Infor         -0.21409732

Para la PC4, VS es la única destacable.

3c: Comentar los gráficos que representan las variables en los planos formados por las componentes, intentando explicar lo que representa cada componente

En la primera visualización se comprobará cuál es la correlación respecto a las 2 primeras componentes:

fviz_pca_var(pca_4,
             axes = c(1,2),
             col.var = "cos2", # Color by contributions to the PC
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )

La mayoría de las variables están muy fuertemente correlacionadas de forma positiva con la componente 1 como son Población, el PIB, Industria, Construcción, etc. Mientras, en la componente 2, la Mortalidad posee una alta correlación negativa con la misma, a diferencia de la Natalidad o la TasaParo que se correlacionan de forma positiva. Por tanto, la componente 1 recoge variables que muestran el número de empresas de diferentes sectores (variables más puramente económicas) a la par que la componente 2 explica variables sociales como Natalidad, Mortalidad, tasas de paro, etc.

En cuanto a las componentes 3 y 4:

fviz_pca_var(pca_4,
             axes = c(3,4),
             col.var = "cos2", # Color by contributions to the PC
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )

En este gráfico hay bastante menos correlación que en el anterior destacando en la componente 3 la representación de indicadores socioeconómicos como el paro, la actividad o el IPC, así como las explotaciones agrarias. En la componente 4 no hay apenas representación aunque podría solaparse con la interpretación de la 3 al formar las flechas casi 45º entre componentes.

3d: Mostrar la tabla y los gráficos que nos muestran la proporción de la varianza de cada variable que es explicado por cada componente. ¿Cuál de las

variables es la que está peor explicada?

Esta pregunta se puede responder comporbando los cos2 de cada variable pues muestra la proporción de la varianza de cada variable que es explicada por cada componente:

knitr::kable(var$cos2, digits = 2, caption="Cosenos al cuadrado")

Cosenos al cuadrado

	Dim.1	Dim.2	Dim.3	Dim.4
Poblacion	0.99	0.00	0.00	0.00
Mortalidad	0.13	0.71	0.06	0.02
Natalidad	0.02	0.63	0.12	0.01
IPC	0.14	0.34	0.11	0.18
NumEmpresas	0.99	0.00	0.00	0.00
Industria	0.94	0.01	0.00	0.00
Construccion	0.99	0.01	0.00	0.00
CTH	0.98	0.00	0.00	0.00
Infor	0.91	0.00	0.01	0.05
AFS	0.98	0.00	0.00	0.01
APT	0.97	0.00	0.00	0.02
TasaActividad	0.15	0.28	0.22	0.20
TasaParo	0.00	0.55	0.25	0.05
Ocupados	0.99	0.00	0.00	0.00
PIB	0.97	0.00	0.00	0.02
CANE	0.00	0.02	0.70	0.07
TVF	0.97	0.00	0.02	0.00
VS	0.34	0.01	0.14	0.30

corrplot(var$cos2, is.corr = FALSE, tl.cex=0.6, tl.col="black",
         cl.ratio = 1)

La componente 1 explica la práctica totalidad de múltiples variables mencionadas en el apartado previo (Industria, Infor, PIB, Ocupados, etc.) y la 2 aquellas no tan explicadas por la 1, de manera que se compensan esos nichos y por ello el poder explicativo es tan alto.

fviz_cos2(pca_4, choice = "var", axes=1:4, tl.cex=0.6)

La variable peor explicada es el IPC, con una varianza total explicada por las 4 componentes de \(0.14+0.34+0.11+0.18 = 0.77\)

3e: Mostrar la tabla y los gráficos que nos muestran la proporción de la varianza de cada componente que es explicado por cada variable ¿Cuál de las

variables es la que está peor explicada?

En este caso se deberá atender al apartado contrib (contribuciones) del objeto var:

knitr::kable(var$contrib, digits = 2, caption="Contribuciones de las variables")

Contribuciones de las variables

	Dim.1	Dim.2	Dim.3	Dim.4
Poblacion	8.62	0.00	0.25	0.29
Mortalidad	1.13	27.79	3.57	2.60
Natalidad	0.17	24.54	7.33	1.20
IPC	1.21	13.35	6.88	18.93
NumEmpresas	8.65	0.07	0.01	0.48
Industria	8.16	0.20	0.22	0.05
Construccion	8.60	0.21	0.01	0.07
CTH	8.58	0.01	0.24	0.08
Infor	7.92	0.18	0.42	4.91
AFS	8.55	0.03	0.16	0.84
APT	8.45	0.09	0.08	2.01
TasaActividad	1.31	10.93	13.16	21.44
TasaParo	0.02	21.30	15.00	4.85
Ocupados	8.67	0.03	0.00	0.36
PIB	8.46	0.13	0.14	1.79
CANE	0.03	0.93	43.13	7.75
TVF	8.50	0.00	1.00	0.19
VS	2.97	0.23	8.42	32.16

corrplot(var$contrib, is.corr = FALSE, tl.cex=0.6, tl.col="black",
         cl.ratio = 1)

Las contribuciones de las variables a la variabilidad de las componentes es algo más escueta que viceversa, destacando fundamentalmente las aportaciones de CANE a la componente 3, mortalidad y natalidad a la 2 y VS a la 4 de forma individual.

fviz_contrib(pca_4, choice = "var", axes=1, tl.cex=0.6)

En la componente 1, ninguna variable explica más del 10% de su variabilidad, siendo la que menos aporta la Tasaparo, confirmando que en esta componente contribuyen más las variables puramente económicas que las socioeconómicas o demográficas.

fviz_contrib(pca_4, choice = "var", axes=2, tl.cex=0.6)

A la componente 2, por su carácter más demógrafico y social, ortogonal a la 1, las que menos aportan son variables como TVF o la Población

fviz_contrib(pca_4, choice = "var", axes=3, tl.cex=0.6)

A la componente 3 contribuyen de forma nula los Ocupados o NumEmpresas

fviz_contrib(pca_4, choice = "var", axes=4, tl.cex=0.6)

En la componente 4 las variables que menos aportan con Construcción e Industria

3f: Sobre los gráficos que representan las observaciones en los nuevos ejes y el gráfico Biplot., teniendo en cuenta la posición de las provincias en el gráfico Comentar las provincias que tienen una posición más destacada en cada componente, en positivo o negativo, ¿Qué significa esto en términos

socioeconómicos para estas provincias?

Para responder a ello se presentarán una serie de biplots con los individuos (las provincias) sobre los ejes de las componentes de la 1 a la 4:

fviz_pca_ind(pca_4, col.ind = "cos2", 
             axes = c(1,2),
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             col.cex=0.2,
             repel = TRUE # Avoid text overlapping (slow if many points)
             ) +xlim(-10,20)

En este gráfico se observa que en la PC1, Barcelona y Madrid tienen un valor muy superior a todas las demás provincias por ser los motores económicos y destacar en altos índices de empresas de todo tipo y población. Todas las demás se encuentran muy cerca del origen de coordenadas por lo que ninguna es excesivamente baja en estos aspectos económicos salvo aquellas más rurales, que conforman la España vaciada como Soria, Teruel, Albacete, etc.

En la PC2, la componente de métricas más sociales, destacan com valores positivos y altos Ceuta y Melilla, Cádiz, Almería, Málaga… en general provincias más del sur que presentan unas tasas de natalidad y paro y superiores a zonas del centro y norte como Zamora, Ourense, Lugo, Cantabria, etc. que poseen una mayor mortalidad por una población más envejecida, pero con menor nivel de desempleo.

Para las componentes 3 y 4:

fviz_pca_ind(pca_4, col.ind = "cos2", 
             axes = c(3,4),
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             col.cex=0.2,
             repel = TRUE # Avoid text overlapping (slow if many points)
             ) +xlim(-10,20)

Madrid y Barcelona ya no se encuentran completamente separadas en la componente horizontal (en este caso la 3), ni tampoco ninguna otra provincia se aleja del origen de coordenadas, mostrando poca relación con las explotaciones agrarias o una tasa de paro fuera de lo normal. En cuanto a la componente 4, ciudades como Madrid, Melilla y Ceuta se encuentran a parte por una mayor natalidad o tasa de paro.

A continuación, se representarán conjuntamente las observaciones con las variables dispuestas en los ejes:

fviz_pca_biplot(pca_4, repel=TRUE, col.var="#2E9FD9",
                col.cex=0.2, col.ind = "#696969")+xlim(-10,20)

En esta representación se puede ver provincias como Lugo, Zamora, Ourense o Teruel tienen una alta mortalidad. Melilla, Ceuta y Cádiz destacan por su natalidad. Sevilla, Murcia y Ávila tienen buena tasa de actividad y las ya mencionadas Madrid, Barcelona y en menor medida Valencia, Alicante y Navarra son los motores económicos.

fviz_pca_biplot(pca_4, repel=TRUE, axes=c(3,4), col.var="#2E9FD9",
                col.cex=0.2, col.ind = "#696969")+xlim(-10,20)

3g: Si tuviéramos que construir un índice que valore de forma conjunta el

desarrollo económico de una provincia, como se podría construir utilizando una combinación lineal de todas las variables. ¿Cuál sería el valor de dicho índice en Madrid? ¿Cual sería su valor en Melilla?

En los anteriores apartados se ha comprobado que la componente 1 es la que reúne las variables económicas, por lo que esta podría ser un buen índice resultado de la combinación lineal de las variables de forma estandarizada. Es decir, la PC1 viene dada por la expresión: $ PC1 = 0,294Poblacion - 0,106Mortalidad + 0,041Natalidad + … + 0,292TVF + 0,172VS $ en la que cada para cada provincia se tomará el valor que esta presenta en alguna de las variables y se le restará la media para el valor de dicha variable, todo ello dividido entre la desviación típica de la variable. En el caso de Madrid por ejemplo hay una población de 6454440 a la que se restará la media de 899448 y se dividirá por 1154297:

mean(data$Poblacion, na.rm = TRUE)

## [1] 899448.9

sd(data$Poblacion, na.rm = TRUE)

## [1] 1154297

Después se multiplicará el resultado por 0,294. Esto se repetiría para cada variable y cada provincia, pero se puede extraer directamente de:

ind<-get_pca_ind(pca_4)
knitr::kable(ind$coord, digits =3,caption = "Valores de los individuo
s en las Cp")

Valores de los individuo s en las Cp

	Dim.1	Dim.2	Dim.3	Dim.4
Albacete	-1.410	0.473	0.096	-0.458
Alicante	3.384	0.540	1.919	2.420
Almería	-0.617	2.614	0.208	-0.142
Álava	-1.444	-0.001	-2.032	-0.113
Asturias	-0.204	-1.953	1.298	-0.714
Badajoz	-1.048	1.168	1.796	-0.854
Balears	1.526	0.260	-2.519	2.364
Barcelona	13.683	-1.612	-0.867	-0.279
Bizkaia	0.576	-1.508	-1.180	-0.383
Burgos	-1.202	-1.550	-0.892	0.989
Cantabria	-0.849	-1.126	-0.594	0.401
Castellón	-0.690	0.821	0.518	0.382
Ceuta	-2.125	3.326	-1.811	-1.165
Ciudad Real	-1.392	0.815	1.819	-0.774
Coruña	0.635	-1.442	0.759	0.244
Cuenca	-2.132	-0.569	1.048	-0.932
Cáceres	-1.503	-0.180	1.269	-0.223
Cádiz	0.128	1.771	0.369	-0.249
Córdoba	-0.594	0.811	1.292	-0.169
Gipuzkoa	-0.281	-1.485	-1.879	0.266
Girona	0.388	0.544	-1.387	1.638
Granada	-0.072	1.124	1.511	0.500
Guadalajara	-1.408	1.542	-1.849	0.906
Huelva	-1.265	1.168	0.013	-0.327
Huesca	-1.776	-0.984	-0.587	0.304
Jaén	-1.221	1.424	3.407	-0.479
León	-1.464	-2.016	0.877	-0.809
Lleida	-0.835	-0.136	-1.472	1.185
Lugo	-1.826	-2.785	0.871	-0.353
Madrid	16.778	-0.366	-0.849	-2.365
Melilla	-2.218	4.782	-1.905	-2.231
Murcia	1.522	1.442	0.580	0.901
Málaga	2.006	1.325	0.869	1.104
Navarra	-0.653	0.078	-1.096	-0.081
Ourense	-1.965	-2.858	1.098	-1.204
Palencia	-2.122	-1.951	-0.695	-0.142
Palmas	0.092	1.857	-0.330	-0.824
Pontevedra	0.036	-0.607	0.052	-0.328
Rioja	-1.383	-0.484	-1.354	0.288
Salamanca	-1.612	-1.425	-0.121	-0.086
Santa Cruz	-0.029	1.573	-0.126	-0.432
Segovia	-1.931	-1.015	-1.110	0.263
Sevilla	1.948	1.775	0.712	-0.546
Soria	-2.399	-1.857	-0.778	-0.503
Tarragona	0.175	1.040	-0.102	1.512
Teruel	-2.185	-1.082	-0.351	-0.413
Toledo	-0.461	1.449	0.812	0.463
Valencia	4.770	0.360	2.961	2.349
Valladolid	-1.007	-0.704	-1.052	0.196
Zamora	-2.287	-3.169	0.578	-0.622
Zaragoza	0.115	-0.237	-0.156	0.068
Ávila	-2.152	-0.982	0.365	-0.545

Madrid tendría un valor respecto a potencia económica de 16,778 mientras que Melilla presentaría un -2,218 haciendo notar numéricamente la diferencia abismal entre ambas que quedó manifiesta en los gráficos.

Pregunta 4: Representar un mapa de calor de la matriz de datos, estandarizado y sin estandarizar para ver si se detectan inicialmente grupos de provincias.

Se comienza creando una nueva versión de los datos (data_st) estandarizados a partir del dataset original y eliminando las variables no numéricas:

data_st<-scale(data)

A continuación se crean mapas de calor para ambos conjuntos:

heatmaply(data, seriate="mean",
          row_dend_left = TRUE, plot_method = "plotly")

heatmaply(data_st, seriate="mean",
          row_dend_left = TRUE, plot_method = "plotly")

Con los datos no estandarizados, parece muy clara la división en 3 grupos del conjunto de datos, Madrid y barcelona por un lado, Las potencias turísticas e industriales del mediterráneo y cantábrico por otro y todas las demás en un tercer grupo (aunque podría hacerse una cuarta segregación dando lugar a 4 grupos). No obstante, al no estandarizar y permitir que las variables con distintas unidades de medida se mezclen (población total y tasa de actividad en % por ejemplo) los resultados no reflejan bien las relaciones.

En el segundo gráfico con estas ya estandarizadas aparece una distinción entre grupos menos clara, lo que puede lleva a deducir que hay más grupos: Madrid y Barcelona, las ciudades autónomas, Valencia y Alicante y luegos 2 grupos adicionales con el resto de comunidades.

Pregunta 5: Realizar un análisis Jerárquico de clusters para determinar si existen grupos de provincias con comportamiento similar.

Para llevar a cabo el análisis cluster, una vez que se han explorado superficialmente los datos con el mapa de calor, se debe crear la matriz de distancias entre las variables (únicamente sobre ellas estandarizadas pues es la manera correcta de hacerlo):

dist_st <- dist(data_st, method = "euclidean")

Con el anterior comando se extrae la distancia euclídea.

A continuación se procede a su visualización:

fviz_dist(dist_st)

También se puede crear la visualización mediante el heatmap visto previamente y agrupando los más parecidos.

heatmaply(as.matrix(dist_st), seriate="mean",
          row_dend_left = TRUE, plot_method = "plotly")

Las provincias más agrupadas según se refleja en estos gráficos son Madrid y Barcelona, muy cerca del origen de coordenadas, luego se distingue un gran grupo que incuye provincias como Huesca, Teruel, Zamora, Lugo, Bizkaia, Alava, etc. Un pequeño grupo con Guadalajara, Girona, Baleares y Tarragona. Ceuta y Melilla en su grupo particular y luego algunas otras.

Por último, se procede a representar el dendrograma de los datos estandarizados mediante el método ward.D2:

res_hc <- hclust(dist_st, method="ward.D2")
fviz_dend(res_hc, cex = 0.5)

5a: A la vista del dendrograma ¿Cuántos clusters recomendarías?

Según el dendrograma, podría haber distinto número de clusters. La opción más adecuada sin sobreestimar es 3 clusters, a la altura de entre 10 y 15, o entrando en cada uno de los 2 grupos más a la derecha podrían darse hasta 6 como se muestra en las siguientes representaciones:

# 3 clusters
as.dendrogram(res_hc) %>% set("branches_k_color", 
             value = c("red", "blue", "green"), k = 3) %>%
    plot(main= "3 clusters")

# 6 clusters
as.dendrogram(res_hc) %>% set("branches_k_color", 
             value = c("red", "blue", "green", "brown", "orange","purple"), k = 6) %>%
    plot(main= "3 clusters")

Se deberá comprobar con los métodos de optimización cuál es el mejor número.

5b: Representar los individuos agrupados según el número de clusters elegido.

La opción más razonable es hacer 3 clusters por lo que se tomará este como número óptimo por el momento. Se comprueba cuántos individuos hay por cluster:

group<-cutree(res_hc, k=3)
knitr::kable(table(group), caption="Nº individuos por cluster")

Nº individuos por cluster

group	Freq
1	31
2	19
3	2

Se visualizan estos grupos representados en el dendrograma:

fviz_dend(res_hc, k=3, cex=0.5, color_labels_by_k = TRUE, rect=TRUE)

También se pueden visualizar los cluster sobre las componentes principales:

fviz_cluster(list(data = data_st,cluster=group),
             ellipse.type = "convex",
             repel = TRUE,
             ggtheme = theme_minimal(),
             show.clust.cent=FALSE,
             main = "Representación sobre PC1 y PC2"
             )

5c: ¿Qué número óptimo de clusters nos indican los criterios Silhoutte y de Elbow?

El número óptimo de clusters según el criterio Elbow y de Silhouette sería:

# Elbow method
fviz_nbclust(data_st, kmeans, method = "wss") +
    geom_vline(xintercept = c(3,6), linetype=2)

  labs(subtitle = "Elbow method")

## $subtitle
## [1] "Elbow method"
## 
## attr(,"class")
## [1] "labels"

# Silhouette method
fviz_nbclust(data_st, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")

El criterio de Elbow muestra 2 puntos de inflexión claros en n=3 y n=6, tal como se mencionaba en apartados anteriores. En principio 3 sigue siendo más estable, pero se deberá probar a visualizar con n=6 por si arroja resultados más adecuados.

Con el criterio de Silhouette ocurre algo similar sólo que los picos se encuentran fundamentalmente en 3 y 5.

5d 1: Con el número de clústeres que nos indica Elbow en el apartado anterior realizar un agrupamiento no jerárquico. Representar los clústeres formados en los planos de las Componentes principales. Relacionar la posición de cada clúster en el plano con lo que representa cada componente principal.

El agrupamiento no jerárquico se llevará a cabo con el algoritmo de kmeans y, antes de proseguir, se comprobará con silhouette la calidad con 3 y 6 clusters:

set.seed(123)
km_res_3<-kmeans(data_st, 3)
km_res_6<-kmeans(data_st, 6)

sil_3 <- silhouette(km_res_3$cluster, dist(data_st))
rownames(sil_3) <- rownames(data_st)

sil_6 <- silhouette(km_res_6$cluster, dist(data_st))
rownames(sil_6) <- rownames(data_st)

fviz_silhouette(sil_3)

##   cluster size ave.sil.width
## 1       1    2          0.66
## 2       2   23          0.08
## 3       3   27          0.41

fviz_silhouette(sil_6)

##   cluster size ave.sil.width
## 1       1    2          0.46
## 2       2   16          0.16
## 3       3   13          0.36
## 4       4    2          0.60
## 5       5    4          0.21
## 6       6   15          0.24

Tanto con 3 como con 6 cluster se obtiene un valor de with muy similar, con algunos valores negativos (incorrectamente asignados en ambos). Al ejecutar varias veces (y por tanto variar la selección aleatoria con kmeans) los resultados muestran grupos poco constantes, a veces sin widths negativos, grupos muy reducidos o muy extensos.

Por ello, se decide crear una función que muestre cuál es la media de las widths según 3,5 o 6 clusters tras ejecutar 100 veces el algoritmo:

get_average_width<-function(clusters){
    widths=0
    counter=0
    for (i in c(1:100)) {
        res<-kmeans(data_st, clusters)
        sil <- silhouette(res$cluster, dist(data_st))
        rownames(sil) <- rownames(data_st)
        widths=widths+mean(sil[,3])
        counter=counter+1
    }
    
    average=widths/counter
    return(average)
}

get_average_width(3)

## [1] 0.2819402

get_average_width(5)

## [1] 0.2297499

get_average_width(6)

## [1] 0.2301892

Tras esta comprobación adicional, 3 clusters ofrecen la mayor width media, de manera que se tomará este número como óptimo para continuar con el análisis.

Se procede a representarlos sobre las componentes principales:

fviz_cluster(list(data = data_st,cluster=group),
             ellipse.type = "convex",
             repel = TRUE,
             ggtheme = theme_minimal(),
             show.clust.cent=FALSE,
             main = "Representación sobre PC1 y PC2"
             ) +xlim(-10,20)

La interpretación de esta representación es la misma que se hacía con el análisis de las componentes. Cuanto más a la derecha en el gráfico, mayor es la concentración empresarial, poblacional y en general, mejor situación económica presentan las provincias. Por otro lado,cuanto más hacia arriba se encuentren las provincias, gozarán de mejor situación social en cuanto a natalidad o tasa de actividad mientras que hacia abajo supone mayor mortalidad.

Sobre PC3 y PC3:

fviz_cluster(list(data = data_st,cluster=group),
             axes = c(3,4),
             ellipse.type = "convex",
             repel = TRUE,
             ggtheme = theme_minimal(),
             show.clust.cent=FALSE,
             main = "Representación sobre PC3 y PC4"
             ) +xlim(-10,20)

Las componentes PC3 y PC4 eran algo menos específicas que las 2 primeras destacando magnitudes como las explotaciones agrarias, la tasa de actividad y el IPC. Esta agrupación de conjuntos muestra que cuanto más a la derecha en el gráfico, más agricultura hay en las provincias como es en los casos de Alicante, Valencia o Murcia y cuanto más a la izquierda y en la parte positiva del eje Y se ubican empresas con buena tasa de actividad y un IPC moderadamente elevado como Girona, Tarragona o Huelva.

5d 2: Evaluar la calidad de los clusters.

Parte de este apartado se realizaó al comienzo del apartado, pero de manera algo más general para averiguar el número de clusters correcto. A continuación se llevará a cabo de forma más minuciosa sabiendo que se dispone de 3 clusters:

set.seed(123)
sil_3 <- silhouette(km_res_3$cluster, dist(data_st))
rownames(sil_3) <- rownames(data_st)
head(sil_3[,1:3])

##          cluster neighbor  sil_width
## Albacete       3        2 0.23404690
## Alicante       2        3 0.12861643
## Almería        2        3 0.23985301
## Álava          3        2 0.33700472
## Asturias       3        2 0.38504449
## Badajoz        2        3 0.03338686

fviz_silhouette(sil_3)

##   cluster size ave.sil.width
## 1       1    2          0.66
## 2       2   23          0.08
## 3       3   27          0.41

Con el primer head se muestra el valor del width que posee cada una de las provincias siendo que cuanto más cerca de 1 esté, mejor agrupado estará el individuo y si está próximo a 0, podría ubicarse en cualquier otro grupo. Con la función para representar el gráfico se muestra la media de esta medida para cada uno de los clusters con sus correspondientes observaciones. Los cluster de 27 y 2 individuos tienen un valor muy elevado, pero el de 23 miembros se ve mermado por las provincias con valor negativo (Girona, Huelva, Guadalajara, Baleares y Castellón y Ciudad Real) que, pese a todo, están muy cerca de 0.

sil_3[,1:3] %>% as.data.frame() %>% filter(sil_width<0)

##             cluster neighbor   sil_width
## Balears           2        3 -0.03186983
## Castellón         2        3 -0.06554316
## Ciudad Real       2        3 -0.05790974
## Girona            2        3 -0.04903299
## Guadalajara       2        3 -0.06677319
## Huelva            2        3 -0.06098097

5e: Explicar las provincias que forman cada uno de los clústeres y comentar cuales son las características socioeconómicas que las hacen pertenecer a dicho clúster.

sorted<-sort(km_res_3$cluster)
knitr::kable(sorted, caption="Provincia y cluster")

Provincia y cluster

	x
Barcelona	1
Madrid	1
Alicante	2
Almería	2
Badajoz	2
Balears	2
Castellón	2
Ceuta	2
Ciudad Real	2
Cádiz	2
Córdoba	2
Girona	2
Granada	2
Guadalajara	2
Huelva	2
Jaén	2
Melilla	2
Murcia	2
Málaga	2
Palmas	2
Santa Cruz	2
Sevilla	2
Tarragona	2
Toledo	2
Valencia	2
Albacete	3
Álava	3
Asturias	3
Bizkaia	3
Burgos	3
Cantabria	3
Coruña	3
Cuenca	3
Cáceres	3
Gipuzkoa	3
Huesca	3
León	3
Lleida	3
Lugo	3
Navarra	3
Ourense	3
Palencia	3
Pontevedra	3
Rioja	3
Salamanca	3
Segovia	3
Soria	3
Teruel	3
Valladolid	3
Zamora	3
Zaragoza	3
Ávila	3

También convendría representar algunos boxplot con los valores que poseen las provincias que conforman los clusters en las variables más correlacionadas con las componentes 1 y 2 (al haber tantas en la variable 1 se elegirán Población y PIB) y para la 2, Natalidad y Mortalidad.

data$Grupo<-as.factor(group)

#Boxplot Población
pob<-ggplot(data, aes(x=Grupo, y=Poblacion, fill=Grupo))+
    geom_boxplot() +
    ylab("")+
    ggtitle("Población")

#Boxplot PIB
pib<-ggplot(data, aes(x=Grupo, y=PIB, fill=Grupo))+
    geom_boxplot() +
    ylab("")+
    ggtitle("PIB")

#Boxplot Natalidad
nat<-ggplot(data, aes(x=Grupo, y=Natalidad, fill=Grupo))+
    geom_boxplot() +
    ylab("")+
    ggtitle("Natalidad")

#Boxplot Mortalidad
mor<-ggplot(data, aes(x=Grupo, y=Mortalidad, fill=Grupo))+
    geom_boxplot() +
    ylab("")+
    ggtitle("Mortalidad")


ggarrange(pob, pib, nat, mor + rremove("x.text"),
          ncol = 2, nrow = 2)

Con estas visualizaciones se confirman las deducciones realizadas a lo largo de todo el análisis. El grupo 3 que abarca a Madrid y Barcelonaes el cluster con mejor situación socioeconómica, con valores mucho más elevados que el resto en cuanto a PIB, Población y otros indicadores económicos, a parte de tener una buena Natalidad (por haber más población).

Pasando al grupo 2, este incluye provincias no tan excepcionales en términos económmicos como los núcleos del país, pero sí por encima de las del grupo 3. Fundamentalmente son provincias del sur y costa mediterránea junto con muy pocas del interior. Estarían en esta posición por el atractivo turístico ,presencia de buenos puertos para el comercio o proximidad con Madrid Y Barcelona como el caso de Baleares, Málaga, Tarragona, Castellón, Toledo, etc.

Finalmente las provincias del grupo 1 pertenecen a zonas de Castilla La Mancha y Castilla León junto con el norte peninsular, destacando por pobres valores en cuanto a Natalidad, alta mortalidad (por población envejecida y el éxodo rural) así como menor población y PIB debido a la menor presencia industrial y turística.